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NUMERICAL  ANALYSIS  IN  FRACTURE  MECHANICS 


A.  S.  Kobayashl 

Department  of  Mechanical  Engineering*  University  of  Washington 
Seattle*  Washington*  USA  98195 


ABSTRACT 

'"'O)  Recent  developments  In  four  numerical  techniques  In  structural  mechanics* 
which  are  used  to  extract  fracture  parameters  for  linear  elastic*  nonlinear 
and  dynamic  fracture  mechanics*  are  reviewed.  Primary  emphasis  Is  placed  on 
the  finite  element  methods  for  determining  two-  and  three-dimensional  (2-D  and 
3-D)  stress  Intensity  factors  In  linear  elastic  fracture  mechanics.  Crack 
opening  displacements  (COD)  and  J-Integrals  for  2-D,  stable  crack  growth* 
ductile  fracture,  and  use  of  elastic  finite  element  method  In  Its  generation 
mode  for  obtaining  dynamic  elastic  fracture  parameters  are  discussed.  The 
second  topic  Is  the  finite  difference  method  for  analyzing  the  elasto-dynamlc 
and  elastic-plastic  dynamic  states  In  fracturing  2-  and  3-D  problems.  The  use 
of  a  super  finite  difference  code  to  study  dynamic  ductile  fracture  using  the 
void  growth  and  coalescence  model  Is  discussed.  The  third  topic  Is  the 
boundary  element  method  which  has  evolved  Into  a  practical  tool  for  numerical 
analysis  In  3-D  linear  elastic  fracture  mechanics.  The  final  topic  Is  the 
updated  alternating  technique*  which  was  merged  with  a  3-D  finite  element  code 
and  together  with  a  break-through  In  Its  analytical  formulation,  has  become  a 
cost-effective  numerical  technique  in  solving  part  and  complete  elliptical 
crack  problems  In  3-D  linear  elastic  fracture  mechanics.  Comparisons  between 
the  J-Integral  of  a  3-polnt  bend  specimen,  the  stress  Intensity  factor  for  a 
surface  flaw  specimen  and  the  dynamic  stress  Intensity  factor  of  a  fracturing 
dynamic  tear  test  specimen  obtained  by  various  Investigators  are  made. 

INTRODUCTION 

Successful  applications  of  linear  elastic  fracture  mechanics  (LEFM)  In 
numerous  postmortem  analyses  of  failures  In  aerospace  structures  of  the  early 
60’s  and  Its  expanded  role  In  design  synthesis  In  the  70’s  required  precise 
knowledge  of  the  stress  Intensity  factors  associated  with  cracks.  Such  stress 
Intensity  factors  of  two-dimensional  (2-D)  and  three-dimensional  (3-D)  cracks* 
which  are  subjected  to  complex  loading  conditions  can  only  be  obtained  through 
numerical  techniques.  As  a  result*  numerical  analysis  of  fracture  mechanics 
problems  became  the  most  active  branche  of  structural  mechanics  In  the  1970’s 
C13.  One  of  the  first  numerical  solutions  In  fracture  problems,  however*  Is 
the  finite  difference  elastic-plastic  result  of  Jacobs  In  1950  C2]  which  was 
followed  by  others  In  the  1960’s  [3  -  5].  Swedlow  et  al  C6]»  on  the  other 


hand*  used  finite  el  went  Method  to  study  a  slellar  problw.  In  essence* 
these  numerical  results  on  the  el asto-plastlc  states  In  cracked  plates  were 
ahead  of  their  tlee  since  no  plausible  ductile  fracture  criterion  was 
available  In  the  early  1960's. 

During  this  period  of  soul  searching*  other  specialized  numerical  tech¬ 
niques*  which  were  specially  designed  to  compute  node  I  (opening  eode)  SIF  In 
2-D  crack  problems*  emerged.  Such  techniques  Include  extensive  studies  using 
the  boundary  collocation  method  t7]»  series  expansion  of  a  complex  mapping 
function  [8]  and  the  method  of  Laurent  series  expansion  [9].  The  numerical 
solutions*  which  were  generated  by  these  special  techniques*  are  still  valid 
today  and  are  listed  In  stress  Intensity  factor  handbooks  and  are  liberally 
quoted  In  literature  as  well.  These  techniques*  howwer*  failed  to  generate 
subsequent  supporters  and  thus  will  not  be  discussed  In  this  paper. 

In  contrast*  the  overwhelming  popularity  of  finite  el  went  method  to¬ 
gether  with  the  growing  acceptance  of  linear  elastic  fracture  mechanics  In 
structural  mechanics  In  the  late  1960's  provided  the  Impetus  for  an  orderly 
development  In  the  use  of  finite  el  went  method  for  determining  stress  Inten¬ 
sity  factors  for  2-D  linear  elastic  fracture  mechanics  CIO].  The  explosive 
developments  In  finite  elwent  method  approaches  to  linear  elastic  fracture 
mechanics  and  also  to  nonlinear  as  well  as  dynamic  fracture  mechanics  of  the 
1970 's  are  documented  In  several  review  papers  and  special  conference  pro¬ 
ceedings  Cl*  11  -  16].  Limited  reviews  of  available  linear  elastic  fracture 
mechanics  computer  software  for  fracture  mechanics  are  given  In  References 
C17,  18]. 

Review  papers  covering  the  other  three  topics  of  finite  difference  meth¬ 
od*  boundary  elwent  method  and  alternating  technique  are  few  and  scarce.  The 
boundary  elwent  method*  however*  has  attracted  a  large  core  of  users  and  Its 
applications  to  fracture  mechanics  have  been  presented  at  numerous  confer¬ 
ences. 

The  purpose  of  this  paper  Is  to  review  the  above  four  numerical  tech¬ 
niques  In  fracture  mechanics*  with  particular  wphasls  on  development  of 
finite  elwent  method  following  the  period  covered  In  Reference  Clj. 

FINITE  ELEMENT  METHOD 

The  above  mentioned  popularity  In  the  use  of  finite  elwent  method  In 
every  aspect  of  fracture  mechanics  has  resulted  In  technical  papers  too 
numerous  to  be  Included  In  this  review.  Thorough  reviews  on  the  applications 
of  finite  elwent  techniques  to  2-D  static  and  quasi-static  problws  In 
fracture  mechanics  through  the  1970's  have  appeared  In  References  Cl*  11  - 
15].  The  historical  and  Important  developments  of  this  era  will  not  be 
repeated  as  this  paper  will  concentrate  on  the  3-D  static  and  2-D  nonlinear 
and  dynamic  analyses  which  werged  during  and  after  this  period. 

A.  3-D  Static  LEFM  Singularity  Elwent 

Although  3-D  finite  elwent  method  codes  are  available  commercially*  the 
mode  I  stress  Intensity  factor  for  a  seemingly  simple  surface  flaw  In  a  uni¬ 
axial  tension  plate  requires  Inordinate  amount  of  computer  time.  The  densely 
packed  3-D  constant  strain  quadrllaterlals  along  the  curved  crack  front  [19] 


-2- 


for  proper  modeling  of  the  l/fr* stress  singularity  results  In  an  Inefficient 
use  of  computer  time.  Although  crack  opening  displacements  (COD)  was  used  to 
Improve  the  accuracy  of  mode  I  stress  Intensity  factor*  K_,  the  accuracy  of 
such  brute  force  amputation  remains  In  doubt.  The  3-D  counterpart  [20]  of 
the  virtual  crack  extension  method  [21]  attempted  to  Increase*  by  computing 
the  local  strain  energy  change  for  small  crack  tip  displacements  In 
conventional  displacement  elements*  the  solution  accuracy  without  excessive 
number  of  finite  elements. 

a.  Singularity  Element 

Computational  efficiency  can  also  be  Improved  by  Incorporating  the  l/fr* 
strain  singularity  In  the  displacement  elements.  Raju  and  Newman  [22]  used 
such  singularity  element  and  reduced  the  effect  of  Interelement  displacement 
discontinuity  by  surrounding  the  crack  tip  with  two  layers  of  "square-root" 
elements.  A  series  of  nodal  forces  adjacent  to  the  crack  tip  was  then  used  to 
compute  K,.  The  multitude  of  3-D  problems  analyzed  by  this  procedure  Include 
the  surface  flaw  problems  [23]  In  pressurized  cylinders. 

b.  Collapsed  Quarter-Point  Isoparametric  Element 

The  popularity  [1]  of  Barsoum's  collapsed  quarter-point  Isoparametric 
element  [24]  Is  due  to  Its  simplicity  In  execution  which  does  not  require 
special  subroutines  to  available  3-D  finite  element  method  codes.  These 
elements  have  the  correct  l//r  singularity  and  together  with  proper  stress 
Intensity  factor  extraction  procedure*  will  yield  stress  Intensity  factor  of 
sufficient  accuracy  along  the  crack  front.  While  many  use  the  crack  opening 
displacement  procedures  or  the  crack-tip  stress  formula  to  compute  the  stress 
Intensity  factor  In  3-D  problems*  few  procedures  are  developed  specifically 
for  3-D  applications.  In  the  following*  two  such  procedures  are  described. 

Ingraffea  [26]  has  shown  that  for  collapsed  20-node  Isoparametric  ele¬ 
ments  surrounding  the  crack  front  of  Figure  1*  the  three  modes  of  stress 
Intensity  factors  can  be  written  In  terms  of  the  mapped  curvilinear  coordin¬ 
ates  of  n  and  £  as 
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Fig.  (1)  -  Quarter-point  Wedge  Elements 
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where  E  and  v  are  the  modulus  of  elasticity  and  Poisson's  ratio*  respective¬ 
ly.  The  three  stress  Intensity  factors  at  n  distance  along  the  crack  front  can 
be  computed  by  the  crack  surface  displacements  of  u*  v  and  w  on  the  collapsed 
element  face  of  £  *  -1.  In  terms  of  the  physical  distance  along  the  z 
coordinate 


n  *  -  (  2z/L2  +  1  )  14) 

where  L,  Is  the  element  length  along  the  crack  front  as  shown  In  Figure  1. 
The  quadratic  variation  In  displacement  field  In  the  original  20-node  element 
also  results  In  a  quadratic  variation  In  stress  Intensity  factor  along  each 
crack  front  element. 

The  accuracy  of  the  stress  Intensity  factor  extraction  procedure  Is  de¬ 
monstrated  by  analyzing  the  embedded  elliptical  crack  problem  solved  analy¬ 
tically  by  Green  and  Sneddon  [27].  Because  of  symmetry,  one-eighth  of  an 
elliptical  crack  In  a  cube  composed  of  23  elements  and  141  nodes  was  analyzed 
using  SAP  IV  [28].  The  maximum  errors  In  the  numerical  results*  when  compared 
against  the  theoretical  results  [27],  were  5*  and  7%  for  a  circular  crack  and 
an  elliptical  crack  with  an  aspect  ratio  of  1.5,  respectively.  At  the  tip  of 
the  minor  axis  where  the  SIF  Is  maximum  In  the  elliptical  crack,  the  error  was 
2X. 


c.  Finite  Element  Hybrid  Method 

Since  a  review  on  the  use  of  finite  element  hybrid  method  In  fracture 
mechanics  appeared  In  Reference  [1],  substantial  progress  has  been  made  In 
Improving  Its  computational  efficiency.  Hybrid  formulation  Is  now  restricted 
to  finite  elements  surrounding  the  crack  front  with  the  crack-tip  singularity 
being  preserved  through  assumed  l//r  stress  or  /r  displacement  field.  Compu¬ 
tational  efficiency  Is  achieved  by  the  general  purpose  finite  element  code 
which  models  the  bulk  of  the  boundary  value  problem.  The  effect  of  Inter-ele¬ 
ment  discontinuity  between  the  singular  and  conventional  elements  Is  minimiz¬ 
ed  through  the  hybrid  formulation.  For  2-D  problems,  the  hybrid  method  can  be 
formulated  such  that  one  element,  which  has  the  assumed  1/vr  singularity  and 
also  an  assumed  compatible  boundary  displacements,  completely  encompasses  the 
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crack  tip  [29].  For  3-0  problems*  the  crack  front  becomes  one  of  the  boundar¬ 
ies  of  the  several  hybrid  elements  which  surround  the  crack  tip.  The  assumed 
stress  hybrid  method  Is  based  on  an  assumed  equlllbrlatlng  stress  field*  which 
contains  the  proper  crack-tip  singularity*  and  on  Independently  assumed  bound¬ 
ary  displacements  [30].  In  terms  of  the  local  coordinates  of  n*  z  and  t»  the 
singular  stress  and  corresponding  displacements  fields  are  represented  as 
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where  r  and  e  are  the  local  polar  coordinates  In  the  plane  normal  to  the 
crack  front.  While  2-D  finite  element  method  codes  can  Incorporate  higher 
order  terms  In  r*  such  as  the  constant*  /r*  r  ...  terms  [31]*  corresponding 
series  forms  of  equations  (3)  and  (4)  are  not  available.  The  assumed  stress 
In  the  crack  tip  hybrid  element  Is  thus  represented  as 
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where  the  first  term  of  the  above  Is  equation  (4)  and  the  second  term  Is  a 


regular  polynomial  of  x»  y  and  z.  For  a  traction  free  crack*  this  second  term 
must  satisfy  the  homogeneous  equilibrium  equations  and  the  traction  free 
boundary  conditions  on  the  crack  surface.  The  assumed  element  boundary 
displacements  contain  equation  (7)  for  all  element  boundaries  which  Intersect 
with  the  crack  front  and  satisfy  Inter-element  compatibility.  When  these 
hybrid  elements  are  merged  Into  a  general  purpose  3-0  code*  equation  (7)  will 
Introduce  as  unknowns*  the  three  stress  Intensity  factors  of  K  (m-  1*2  or  3) 
In  addition  to  the  unknown  generalized  nodal  displacement  of  q  as 
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q  and  K  are  thus  obtained  by  solving  equation  (8). 

■\i  m 

Plan  and  Morlya  [32]  used  a  twelve  node*  assumed  stress  hybrid*  half 
element*  as  shown  In  Flnure  2.  to  determine  the  distribution  of  stress 


8-  NODE  BASIC  12- NODE  BASIC  16-NODE  BASIC 
ELEMENTS  ELEMENTS  ELEMENTS 


12-NODE  HALF  20-N0DE  HALF  26-NODE  HALF 
ELEMENT  ELEMENT  ELEMENT 


20-N0DE  36-NODE  46 -NODE 

SUPERELEMENT  SUPERELEMENT  SUPERELEMENT 

Fig.  (2)  Three-Dimensional  Hybrid  Elements 

Intensity  factor  along  a  semi-circular  surface  crack  In  a  tension  specimen. 
Because  of  four  way  symmetry*  only  one  quadrant  of  the  specimen  was  analyzed 
with  284  nodes  and  652  degrees  of  freedom.  The  numerical  results*  when 
compared  with  Smith’s  results  [33]  obtained  by  the  alternating  technique*  were 
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within  12 %  at  the  free  surface  and  coincided  at  the  maximum  crack  depth. 

The  assumed  displacement  hybrid  method  Is  based  on  an  assumed  displace¬ 
ment  field*  which  Includes  equation  (6)  with  K  »  and  Independently  assumed 
element  boundary  displacements  and  boundary  tra&lons.  The  general  format  of 
the  final  element  equation  Is  Identical  to  equation  (8)  where  K  Is  treated 
also  as  three  unknowns  along  the  the  crack  front.  m 

Atlurl  and  Katherlsean  [34]  used  static  condensation  to  produce  a  20  node 
quadratic  Isoparametric  super-element*  as  shown  In  Figure  3*  based  on  the 


Fig.  (3)  Super  Singular  Element 

assumed  displacement  hybrid  method.  Using  this  super-element  In  a  general  3-D 
finite  element  code*  they  analyzed  the  surface  flaw  problem  with  280  finite 
elements  and  4815  degrees  of  freedom.  The  results  are  discussed  In  the  last 
section  entitled  "Benchmark  Solutions"  of  this  paper. 

B.  2-D  Nonlinear  Singularity  Element 

Since  the  early  use  [6]  of  constant  strain  elements  for  elastic-plastic 
finite  element  method  analysis  of  2-D  fracture  problems*  this  simple  element 
Is  still  used  today  with  success*  primarily  due  to  the  moderate  or  lack  of 
stress  singularity  at  the  plastically  yielded  crack  tip.  Kannlnen  et  al  [35] 
used  such  finite  element  method  code  In  both  Its  generation  and  application 
modes  to  study  stable  crack  growth  and  Instability  of  A533-B  steel  and  2219- 
TB7  aluminum*  center  cracked  and  compact  specimens.  Loadllne  displacement  and 
crack  length  measurements  were  Input  Into  generation-mode  calculations  and  the 
applied  load  among  others  were  output  for  evaluation.  Figure  4  shows  the 
computed  and  measured  applied  load  versus  loadllne  displacement  relation  of 
the  steel  compact  specimen.  Other  uses  of  the  conventional  element  Include 
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Fig.  (4)  Load-displacement  Curve  for  A533B  Steel  CT  Specimens 

the  updated  Lagranglan  finite-deformation,  finite-element  analysis  by 
McMeeklng  and  Parks  [36]  who  studied  the  Influence  of  crack-tip  blunting  on 
the  J -based  characterization  of  the  crack  tip  region. 

For  a  material  with  a  strain  hardening  Index  of  n,  the  strain  field  of 
the  dominant  singularity  at  the  crack  tip  under  deformation  theory  of  plasti¬ 
city  Is  [37,  38] 
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where  J  Is  the  J-Integral  as  per  Rice  [39],  a..  and  e. .  are  functions  of  Band 
n,  and  oQ  Is  the  yield  stress.  J  1J 


Early  uses  of  the  above  slngluar  stress  and  strain  fields  Include  that  of 
Hilton  C40]  who  used  deformation  theory  of  plasticity  to  study  ductile  crack 
Initiation  under  monotonlcally  Increasing  load.  Shlh  [40]#  on  the  other  hand# 
constructed  a  circular  crack  tip  element  with  the  dominant  stress  and  strain 
singularities  and  studied  the  changes  In  elastic- pi  as tic  boundaries  with 
variations  In  n  values  under  small  scale  yielding.  Atlurl  et  al  [42]  Incorp¬ 
orated  the  above  stress  and  strain  singularities  Into  the  hybrid-displacement 
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Fig.  (5)  Translation  of  Singular  Element 

finite  elements  surrounding  the  crack  tip.  The  virtual  work  equation  for  an 
Incremental  crack  growth#  Aa#  as  shown  In  Figure  5#  was  used  to  simulate  crack 
closure  and  opening  under  a  cyclic  loading  condition  with  a  single  overload  as 
shown  In  Figure  6.  The  penalty  function  and  superposltlton  method  used  by 
Yagawa  et  al  [43]  Is  similar  In  formulation  where  the  penalty  coefficients  are 
stiffness  coefficients  which  optimize  the  potential  function  under  nonlinear 
constraints. 

For  an  Ideally  plastic  material#  n  ■  »  In  equat1ons(9)  and  (10)  yield 
°I1  *  °o  and  e1i  *****  an  1/r  singularity.  Barsoum  [44]  showed  that  the  same 
triangular  quarter-point  element  [24]  of  elastic  analysis  will  possess  the  1/r 
strain  singularity  when  the  condensed  crack  tip  nodes  are  allowed  to  slide. 
Shlh  et  al  [45]  used  this  quarter  point  element  to  model  crack  tip  blunting 
and  growth.  As  shown  In  Figure  7#  crack-tip  blunting  Is  modeled  by  separating 
the  condensed  crack  tip  nodes  and  crack  extension  Is  modeled  by  sequential 
shifting  of  the  crack  tip  node.  The  crack  tip  opening  displacement  (CTOA)  and 
crack  tip  displacement  can  be  determined  directly  from  Figure  7(c). 
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C.  3-0  Non-linear  Singularity  Element 

The  singularity  element  used  In  3-D  elastic-plastic  analysis  of  cracks 
apparently  Is  limited  to  Barsoum's  triangular  quarter-point  element.  Although 
this  element  contains  only  a  1/r  singularity*  Benz ley  [46]  showed  heurlstlcal- 
1y  that  when  used  with  a  stress-strain  relation  governed  by  a  power  hardening 
law*  this  element  provided  the  correct  singular  states  of  strain  and  stress 
Interior  of  the  element.  deLorenzI  [47]  used  556  of  such  20-noded  Iso¬ 
parametric  elements  In  an  8300  degree- of- freedom  system  to  study  the  elastic- 
plastic  behavior  of  a  surface  flaw  In  the  beltllne  region  of  a  pressurized 
reactor  vessel.  The  semi-elliptical  crack  front  was  surrounded  with  triangu¬ 
lar  quarter-point  elements  and  a  Ramberg-Osgood  power  hardening  law  for 
stress-strain  relation  was  used.  COD  at  the  symmetry  plane*  as  shown  In 
Figure  8*  demonstrated  the  need  for  a  3-D  elastic-plastic  analysis  since  the 
plane  strain  (2-D)  approximation  clearly  overestimated  the  severity  of  the 
flaw. 


Fig.  (8)  COD  of  a  Surface  Flaw 
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D.  2-0  Dynamic  Singularity  Element 

Dynamic  fracture  analysis  has  also  been  conducted  successfully  with 
conventional  elements  during  the  past  five  years.  Such  analysis  Include  the 
extensive  code  verifications  by  Kobayashl  et  al  [48*  493  and  the  recent  work 
by  Jung  et  al  [503.  Accuracy  of  such  simplified  finite  element  method  codes 
can  be  Improved  by  using  a  proper  crack  node  release  mechanism  which  was  the 
subject  of  considerable  debate  In  the  late  1970's  [51  -  533.  This  writer's 
experiences*  however*  Indicate  that  the  simple  node  release  mechanism  of  lin¬ 
early  decreasing  crack- tip  nodal  force  yield  dynamic  stress  Intensity  factor 
which  are  In  good  agreement  with  those  obtained  by  photoelasticity  [543. 
Yagawa  et  al  [553*  on  the  otherhand*  represented  the  crack  surface  traction* 
acting  on  the  extending  crack  suface*  as  a  Lagrange  multiplier*  and  optimized 
the  dissipated  surface  energy  during  crack  extension. 

While  singularity  elements  In  dynamic  finite  element  method  was  used 
earlier  by  Anderson  et  al  [563  and  Aokl  et  al  [573*  the  most  successful  use  of 
such  element  Is  by  Atlurl  et  al  who  Incorporated  the  l//r  singularity  In  the 
displacement  hybrid  crack  tip  element  [58*  593.  Figure  9  shows  the  Atlurl 's 

TYPE  A  ••  MOVING  SINGULAR  ELEMENT 
TYPE  B  :  DISTORTING  REGULAR  ELEMENT 
TYPE  C  :  NON -DISTORTING  REGULAR  ELEMENT 
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Fig.  (9)  Dynamic  Crack  Growth  Using  a  Simplicity  Element 

singular  element*  which  moves  with  the  crack  tip  and  retains  Its  shape*  and 
the  continually  distorting  regular  elements  surrounding  the  crack  tip. 
Periodic  mesh  readjustment  Is  necessary  with  crack  extension  when  this 
procedure  Is  used.  Similar  approach  without  the  distorting  regular  elements 


was  used  by  Gunther  et  al  [60].  Figure  10  shows  good  agreement  between  the 
numerically  generated  dynamic  stress  Intensity  factor  by  Atlurl  et  al  and  the 
theoretical  solution  by  Broberg  [61]. 
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Fig.  (10)  Dynamic  Stress  Intensity  Factor  of  Broberg  Crack 


Literature  Is  void  with  3-D  dynamic  finite  element  analysis*  due  In  part 
to  the  enormity  In  computational  requirements*  but  also  due  to  the  lack  of 
definitive  experimental  observations  on  3-D  dynamic  crack  extension.  Such 
dynamic  crack  extension  history  Is  necessary  to  execute  a  dynamic  finite 
element  code  In  Its  generation  mode.  3-D  dynamic  analysis  by  the  propagation 
mode*  on  the  otherhand*  requires  aprlorl  a  3-D  dynamic  crack  propagation  law 
which  Is  equally  missing  at  this  time. 

FINITE  DIFFERENCE  METHOD 

As  mentioned  In  the  Introduction*  finite  difference  method  predates  the 
now  popular  finite  element  codes  In  Its  application  to  fracture  problems. 
Subsequent  development  of  the  finite  difference  codes  were  taken  over  by  weap¬ 
on  researchers  and  re-emerged  as  vastly  superior  general  purpose  codes  which 


could  solve  3-D  problems  ranging  from  dynamic  plasticity  to  gas  dynamics.  For 
example*  the  HEMP  3-D  code  [62]  Is  an  explicit  finite  difference  code  which 
does  not  require  large  computer  storage.  Crack  tip  singularity  Is  thus  han¬ 
dled  by  swamping  the  crack  region  with  large  number  of  zones  and  no  known  at¬ 
tempt  has  been  made  to  Incorporate  crack  tip  singularities  Into  the  computa¬ 
tion.  Such  simplified  model  results  In  enormous  computer  time  which  normally 
cannot  be  executed  outside  of  special  laboratories.  The  code  can  also  handle 
ductile  fracture  with  relative  ease  due  to  the  reduced  severity  In  stress 
singularity. 

A.  Static  Analysis 

Historically*  the  available  general  purpose  finite  difference  codes  were 
designed  to  solve  complex  dynamic  problems  and  thus  no  static  finite  differ¬ 
ence  programs  for  analyzing  static  fracture  problems  exist  to  date.  In  an 
overkill  attempt  to  demonstrate  the  versatility  of  such  supercodes*  Chen  [633 
showed  that  the  static  stress  Intensity  factor  can  be  obtained  from  a  simple 
average  of  the  stress  waves  set  up  In  the  crack-tip  region  when  high  artifi¬ 
cial  viscosity  Is  Inserted  to  damp  out  the  stress  waves  In  the  3-D  supercodes. 
A  similar  converge  scheme  using  dynamic  relaxation  was  used  by  Shmuely  et  al 
[643  to  obtain  the  static  SIF  In  a  finite  thickness*  central  crack  tension 
plate. 
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Fig.  (II)  Dynamic  Stress  Intensity  Factor  In  an 
Impulse  Loaded  Plate 

B.  Dynamic  Analysis 

For  elastic  analysis*  a  simple  extrapolation  scheme  Is  used  to  extract 
stress  Intensity  factors  from  the  numerically  determined  crack  tip  stresses 
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[63].  Figure  11  shows  the  dynamic  stress  Intensity  factor  In  a  center  cracked 
plate  subjected  to  sudden  tension  loading.  Also  shown  for  comparison  Is 
Baker's  solution  [65]  for  a  pressurized  semi-infinite  crack  which  suddenly 
appeared  In  an  Infinite  plate.  Complex  fracture  problems*  such  as  an  Inter¬ 
nally  or  externally  surface  flawed  cylinder  subjected  to  sudden  pressuriza¬ 
tion  has  also  been  solved  [66]  by  the  HEMP  code.  The  advantage  of  such 
supercode*  however*  lies  In  Its  ability  to  analyze  elastic-plastic  dynamic 
fracture  problems*  such  as  a  notched  bar  subjected  to  sudden  tensile  loading 
[66*  67]  and  a  Charpy  V-notched  specimen  [68].  In  these  ductile  fracture 
analyses*  a  void  growth  and  coalescence  criterion  [68]  was  used  to  predict  the 
onset  and  propagation  of  a  ductile  crack. 

Special  purpose  finite  difference  codes  have  been  used  to  analyze  dynamic 
fracture  problems  but  the  earlier  analyses  [70*  71]  did  not  focus  on  viable 
dynamic  fracture  parameters*  prlmarly  due  to  the  undeveloped  state  of  science 
In  dynamic  fracture  at  that  time.  Finite  difference  method  was  used  by  Pope- 
lar  et  al  [72]  to  study  the  dynamic  elastic  response  of  an  Internally  cracked 
cylinder  subjected  to  Impulse  loading.  The  dynamic  stress  Intensity  factor 
was  determined  by  the  energy  release  rate  calculated  from  the  displacements  at 
the  nodes  In  the  vicinity  of  the  propagating  crack  tip.  Shmuely  et  al  [73] 
Incorporated  a  moving  substructure*  as  shown  In  Figure  12  with  the  proper  1/fr 


Fig.  (12)  Stationary  and  Moving  Grid 

singularity  In  their  finite  difference  program.  One  half  reduction  In  comput¬ 
ing  time  In  static  analysis  and  Improved  accuracy  In  energy  release  rate  comp¬ 
utation  In  dynamic  analysis  are  claimed.  As  a  straightforward  application  of 
a  dynamic  finite  difference  shell  code*  Emery  et  al  [74]  studied  the  axial 
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cracking  of  a  pressurized  crack  with  coupled  depressurization  calculation. 
This  analysis  was  extended  to  large  scale  yielding  with  critical  CTOA  as  a 
dynamic  fracture  criterion  [75]. 

BOUNDARY  ELEMENT  METHOD 

While  the  boundary  element  method*  which  was  earlier  referred  to  as  the 
boundary  Integral  equation  method*  as  a  stress  analysis  tool  dates  back  to 
the  mid  1960's  [76]*  Its  application  towards  solving  fracture  mechanics  prob¬ 
lems  was  pioneered  by  Cruse  [77*  78]  In  the  early  1970's.  Boundary  equation 
method  requires  only  the  discretization  of  the  boundary  of  the  structure  In 
contrast  to  the  domain  discretization  required  by  finite  element  method.  Thus 
boundary  equation  method  Is  not  suitable  for  analyzing  elastic-plastic  and 
el asto-dynaml c  fracture  problems  although  some  recent  work  [79*  80]  suggests 
that  such  use  of  boundary  equation  method  may  not  be  far.  On  the  other  hand* 
the  reduced  system  of  equation  makes  It  suitable  for  solving  3-D  linear 
elastic  fracture  mechanics  problems  although  It  yields  a  fully  populated* 
nonsymmetrlc  system  of  equations.  Also*  for  an  extending  crack*  boundary 
equation  method  requires  only  the  recomputatlon  of  nodes  along  the  crack 
surface  while  finite  element  method  requires  complete  remeshing  around  the 
original  and  extended  crack.  Recent  applications  of  boundary  element  method 
to  fracture  mechanics  has  been  reviewed  by  Cruse  [81]. 

A.  2-D  Boundary  Element  Meythod 

Since  boundary  element  method  emerged  as  a  numerical  tool  during  the 
period  of  successful  application  of  2-0  finite  el ement  method  to  elasto-statlc 
problems  In  linear  elastic  fracture  mechanics*  Its  use  In  2-D  linear  elastic 
fracture  mechanics  did  not  grow  until  a  special  Green's  function  approach  [78* 
82]  was  developed.  Utilizing  the  eleastlc  Green's  function  for  an  Infinite 
plate  with  a  crack  and  an  Interior  point  load*  the  stress  Intensity  factors 
can  be  represented  by  the  following  set  of  path  Independent  Integrals 
enclosing  a  crack  tip  as 

<v  *n>  -f8R!’IS)ui(Q)ds  + 1,  L!,II(zQ)ti(Q)d5  (id 


where  u.  (Q)  and  t.  (Q)  refer  to  the  displacements  and -surface  tractions  at 
boundary  point  G.  Details  of  functional*  R.  and  L,1  »  are  given  In  [83]. 

By  taking  advantage  olf  symmetry*  Cruse  used  only  seven  boundary  points  to 
obtain  Kj  for  a  central  crack  fracture  specimen  with  6.4  CPU  seconds. 

Blandford  et  al  [84]  eliminated  the  complex  arithmetic  Involved  In  the 
above  by  using  traction  singular  quarter-point  boundary  elements  along  each 
side  of  the  crack  tip.  The  procedure  Is  the  boundary  element  method  counter¬ 
part  of  Barsoum's  finite  element  method  procedure  [24]  where  the  midpoint  In 


the  Isoparametric  quadratic  boundary  element  Is  shifted  to  the  quarter  point. 
This  totally  numerical  procedure  eliminates  the  non-uniqueness  of  flat  crack 
modeling  [823  and  the  matrix  singularity  problem  [853.  Atkinson  et  al  [863 
used  the  path-independent  F.  and  M  Integrals  [873  to  determine  K...  In  an  edge 
cracked  square  torsion  bar  with  a  quadratlcal ly  varying  shear  moaQTus. 

B.  3-D  Boundary  Element  Method 

Since  the  first  applications  of  boundary  element  method  to  3-D  linear 
elastic  fracture  mechanics  In  the  early  1970's  [85],  significant  Inprovements 
has  been  made  In  the  computational  algorithm.  The  use  of  quadratic  isopara¬ 
metric  boundary  element  elements  with  quarter  point  nodes  [893  provided 
accurate  COD'S  which  were  used  to  extract  the  K.  values.  Recent  solutions  to 
3-D  linear  elastic  fracture  mechanics  Include  the  surface  flaw  solutions  In  a 
tension  plate  [903  and  In  a  pressurized  cylinder  [913  by  Hellot  et  al  and 
crack  growth  studies  of  surface  flaw  by  Cruse  [923. 

ALTERNATING  TECHNIQUE 


Application  of  the  alternating  technique  to  3-D  linear  elastic  fracture 
mechanics  was  first  Introduced  by  Smith  et  al  [933  who  solved  the  semi-circu¬ 
lar  surface  flaw  problem.  The  procedure  was  extended  to  an  elliptical  crack 
by  Shah  et  al  [943  who  determined  K,  of  an  embedded  elliptical  crack  near  a 
free  surface.  Later#  Smith  et  al  [w53  extended  the  solution  procedure  to  K.. 
and  K...  determination  of  an  elliptical  crack.  The  alternating  techniques  of 
thosexaays  were  relatively  Inaccurate  due  to  the  limited  curve  fitting  capa¬ 
bilities  of  the  third  order  polynomial  of  the  elliptical  crack  pressure  and  of 
the  modeling  of  the  surrounding  finite  geometry.  A  major  breakthrough  In  the 
latter  was  made  by  Browning  et  al  [963  and  Kullgren  et  al  [973  who  used  a  3-D 
finite  element  code  to  model  the  surrounding  finite  geometry.  Grandt  [983  and 
Barrachln  et  al  [993  used  this  procedure  to  analyze  surface  flaw  problems  and 
In  particular#  the  well  studied  surface  flaw  at  a  hole. 


Another  significant  Improvement  was  made  with  the  derivation  of  the 
complete  analytical  solution  [1003  for  an  embedded  elliptical  crack  In  an 
Infinite  solid  and  subjected  to  modes  I#  II  and  III  crack  tip  deformation. 
Nlshloka  et  al  used  twelve  terms  of  a  fifth  order  polynomial  of  this  analy¬ 
tical  solution  together  with  a  standard  3-D  finite  element  code  to  analyze 
among  other#  the  surface  flaw  problem  [1013  and  the  Internally  and  externally 
flawed  pressui  Ized  cylinders  [1023.  Figure  13  shows  the  reduction  In  residual 
stress  with  alternating  cycles  of  Iteration  for  analyzing  an  externally  flawed 
pressurized  cylinder.  Figure  14  shows  the  resultant  stress  Intensity  factor 
which  Is  compared  with  those  of  [103#  1043.  The  finite  element  method  portion 
of  this  alternating  technique  used  96  20-noded  Isoparametric  elements  with 
1815  degrees  of  freedom  and  the  total  CPU  time  was  about  1000  seconds  with  a 
CYBER  74. 


RESIDUAL  STRESS 


Fig.  13  Residual  Stress  on 
Crack  Surface 


Fig.  (14)  Stress  Intensity  Factors 
for  an  Outer  Surface  Flaw 
In  a  Pressurized  Cylinder 
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BENCHMARK  SOLUTIONS 

In  most  cases*  numerical  techniques  are  validated  by  comparing  the  numer¬ 
ical  results  with  known  theoretical  solutions.  By  nature*  these  theoretical 
problems  In  fracture  mechanics  are  simple  In  geometry  and  uncomplicated  In 
loading  and  thus  agreement  between  the  numerical  and  theoretical  solutions  are 
bound  to  be  good.  Theoretical  solutions  for  more  realistic  fracture  problems 
do  not  exist  and  a  consensus  between  various  numerical  solutions  does  not 
guarantee  their  correctness.  Lacking  other  means  of  comparative  study*  bench¬ 
mark  solutions  which  are  well  defined  boundary  and  Initial  value  problems*  are 
used  to  eliminate  any  ambiguity  In  problem  definition  and  to  assure  that  all 
numerical  solutions  relate  to  the  same  problem.  Three  such  benchmark  problems 
are  discussed  In  the  following. 

A.  2-D  Elastic-Plastic  Crack  Problem 

In  1975*  ASTM  Committee  E24.01.09  undertook  a  task  to  compare  numerical 
solutions  to  elastic-plastic  plane  strain  problems.  A  three-point  fracture 
toughness  test  specimen  with  a  uniaxial  stress-strain  relation  of  A533B  steel 
was  analyzed  by  10  respondents  and  the  edited  and  assembled  solutions  were 
presented  by  Wilson  [105].  Figure  15  shows  the  average  J-value*  which  was 


computed  over  several  Integration  paths  by  each  Investigator*  with  Increasing 
applied  load.  While  substantial  progress  has  been  made  In  the  elastic-plastic 
codes  In  the  ensuing  six  years*  the  wide  differences  In  the  results  obtained 
In  the  late  1970's  are  still  Indicative  of  the  lack  of  consensus  for  valid  2-D 
elastic-plastic  codes  as  well  as  for  basic  physical  laws*  such  as  the  consti¬ 
tutive  relations  under  plastic  flow  under  severe  strain  gradients. 
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B.  3-D  LEFM  Problam 


The  surface  flaw  problam  Is  an  outgrowth  of  a  1976  workshop  at  Battel 1e 
Colunbus  Laboratories  at  which  ties  several  3-D  benchMrk  problems  were  desig¬ 
nated  for  numerical  analysis.  The  numerical  results  for  the  saml-ell Iptlcal* 
surface  flaw  plates  In  tension  and  bending  were  subsequently  assembled  and 
edited  by  McGowan  [106].  Figure  16  shows  the  normalized  stress  Intensity 
factors  by  six  Investigators  where  moderate  differences  are  seen. 
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Fig.  (16)  Stress  Intensity  Factor  for  a 
Surf  ace- Flawed  Tension  Plate 

C.  2-D  Dynamic  Crack  Problem 

While  no  formal  benchmark  problem  In  2-D  dynamic  fracture  was  ever  desig¬ 
nated  for  round  robin  studies*  an  Informal  comparative  study  was  made  between 
the  Investigators  at  the  Georgia  Institute  of  Technology*  Battel le  Columbus 
Laboratories  and  the  University  of  Washington.  Dynamic  crack  propagation  In 
an  A533B  steel*  dynamic  tear  test  (DTT)  specimen  was  analyzed  by  generation 
calculation  using  the  crack  velocity  data  provided  In  [107].  Figure  17  shows 
the  dynamic  SIF  computed  by  the  three  Investigators.  All  three  elastic 
results  are  In  good  agreonent  with  each  other. 
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Fig.  (17)  Dynamic  Stress  Intensity  Factor  of  a 
Dynamic  Tear  Test  Specimen 

CONCLUDING  REMARKS 

The  concluding  remarks  In  this  paper  could  have  been  taken  from  Reference 
Cl]  since  the  four  cited  and  the  numerous  uncited  numerical  techniques  In 
fracture  mechanics  continue  to  develop  at  equal  or  even  higher  pace  since 
Reference  Cl]  was  written.  Also*  the  efficiency  In  numerical  techniques 
continues  to  improve  despite  the  rapid  Increase  In  the  complexity  of  the 
problems  studied. 

Despite  this  explosive  rate  of  development*  coordinated  efforts  to  evalu¬ 
ate  the  old  as  well  as  new  numerical  codes  through  benchmark  problems  are 
scarce.  In  order  to  properly  assess  the  numerous  codes  available  today*  an 
International  effort  In  code  verification  Is  badly  needed  at  this  time. 
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11  October i,  VA  24061 

Profaaaor  A.  C.  triages 
Princeton  UaLvoralty 

Department  of  Aeroapoco  nod 
Mac  ton  leal  Sclancaa 
Prlocacoo,  KJ  08340 

Profaaaor  ft.  H-  Loo 

Stanford  Uotvaroity 

Division  of  tog Inner tag  Msc tonics 

Stanford ,  CA  34305 

Profasaor  Albert  l.  King 
Wayne  State  Unlveralty 
biomechanics  Raaaarch  Cantor 
Detroit,  Ml  48202 

Dr.  V.  t.  Hodgaon 
Wayne  State  University 
School  of  Medicine 
Detroit,  MI  48202 

Dean  8.  A.  loley 
Nerehwoetorn  University 
Department  of  Civil  Eoglaeorlag 
Evans too,  IL  40201 


Universities  ( Con ' t . ) 

Or.  Sxiutl  ft.  ftatdorf 
University  of  California 
School  of  Engineering 
and  Applied  Science 
Los  Angeles,  CA  90024 

Profeaaor  Isaac  Pried 
Boston  Unlveralty 
Department  of  Mathematics 
Boston,  HA  02213 

Professor  C.  Krvpl 
laaaselaer  Polytechnic  Institute 
Division  of  Engineering 
Englnaarlng  Mac tonic a 
Troy.  Iff  12181 

Or.  Jack  R.  Vtaaom 
University  of  tolawra 
Depart  sent  of  Mechanical  and  An  reaps an 
Englnaarlng  and  the  Cantor  far 
Composite  Materials 
NS  mark,  fit  13711 

8r.  J.  Daffy 
Sreaa  University 
Otvlalaa  af  tag  lamer lag 
Providence,  ftl  08812 

Or.  J.  L.  toad  lorn 
Carnegie -Mel lam  University 
Department  of  Maetonleal  Eaglaaarlag 
Plttabergh,  PA  1321) 
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Universities  (Coa’t.) 

Or.  Jackaom  C.  S.  Tang 
University  of  Maryland 
Da  pa r toon t  of  Mechanical  Engineering 
College  Park,  KD  20742 

Profaaaor  T.  Y.  Chang 
Oalvaralty  of  Akroa 
Dapartmane  of  Civil  laglnoorlag 
Akron,  OH  44)25 

Profeaaor  Charles  w.  Bart 
Unlveralty  af  Oklahoma 
School  of  Aaroapaca,  Mactoalcal, 
and  Hue  1  ear  Ei*  Inter  tag 
tormaa ,  OK  7)013 

Profaaaor  8atya  H.  Atlurl 
Oaorgla  Institute  of  Tachmolagy 
School  of  fagloaarlag  and 
Mocton lea 

Atlanta ,  a  303)2 

Profaaaor  Graham  f.  Corny 
Uolvoratty  af  Tamo  at  Austin 
Dapar tmaat  af  Aaroapaca  IRgl soar tag 
and  lagloaarlng  Mac ton lea 

toot la,  TX  78711 

Or.  8.  8.  Wang 
Oalvaralty  af  lllloala 
Dapar tmaat  of  fhaoratlcal  and 

Applied  Mae  ha*  lee 
Ur be an,  IL  61801 


Or.  V.  K.  Varodaa 

Ohio  state  Unlveralty  ftaoeacsk  Nmadatlem 
Department  of  loglaoortai  toe ban la a 

Colaotoa,  00  4)210 

Or.  I.  tonkin 

Oalvaralty  af  Paaaaylvamla 
hfarmat  af  totaling  mi 
Mat  or  laid  telmnaa 
Callga  af  toglaaart*  ami 
toyl  lad  to  lama 
Philadelphia,  PA  19104 


ladmatrr  and  taaaarch  laatltvtta 


Or.  toman  tohka 
Pamee  AvtDyns 
Otvtoftom  af  lames 
to I ancon  Caifaritlaa 
Bar  ling  tom,  lift  0110) 

Argiaaa  tot tana l  Laboratory 
Ubret?  Sarvlcaa  Dapartmane 
•700  Snath  Case  dvamma 
Argamma,  IL  80440 


Ualvaraitlaa  (Can't.)  , 

ftafmmaar  V.  C.  todga,  Jr. 
Oftivaralty  of  Miamoaota 
Daparcmoat  af  Aaroapaca  O^laaarlag 
ami  Maehamtea 
Mianoa polls.  Ml  354)5 

to.  0.  C.  touekar 
toi varsity  of  Illinois 
toam  of  Easlmoerlqg 
Urbona,  a  61801 

Profaaaor  0.  J.  tomaork 
Oalvaralty  of  Ullmola 
toportmoat  of  Civil  Kaglaamrldg 
Orhamo,  JL  6100) 

Profeaaor  I.  ftalaomar 

Oalvaralty  of  Collforola,  So*  Olmga 

topartasnt  of  Applied  Mechanics 

U  Jolla,  CA  320)7 

Profaaaor  William  4.  took 
Oalvaralty  of  toaoochnaoeto 
toportmoat  of  Maehamtea  and 
Aaroapaca  Bag  laser 

to her at,  MA  01002 

Profaaaor  G.  Herrmann 
Stanford  Uni varsity 
Department  of  Applied  Machaalca 
Stanford ,  CA  94303 

Profeaoor  J.  0.  Actonbsch 
Nor t hues;  University 
toportmoat  of  Civil  Eogloasriog 
Evanston,  IL  60201 

Profaaaor  8.  8*  Dong 
Univaralty  of  California 
toportmoat  of  Machanlea 
Los  Angolan,  CA  90024 

Profaaaor  lure  Paul 
Unlveralty  of  Poanaylvonia 
Tbwne  School  of  Civil  end 
Mechanical  Eng  Inner lag 
Philadelphia,  PA  l 9104 


Industry  and  ftasearch  Institutes  (Coa't 

Dr.  N.  C.  Jungsr 
Cambridge  Acoustical  Associates 
34  ftlndge  Avenue  Extension 
Cambridge,  HA  02140 

Or.  V.  God  1 no 

General  Dynamics  Corporation 
Electric  Soot  Division 
Groton,  CT  06340 

Dr.  J.  E.  Greanspon 

J.  G.  Eng Inner lag  to  anarch  Aaaaclataa 

38)1  Mania  Drive 

Baltimore,  MD  21215 

towport  tow  Shiphuildlag  ami 

toy  OpcIi  Company 
Library 

to mport  tow,  34  2)601 

to.  W.  P.  ftoalch 
McDonnell  Douglas  Corporation 
5301  Bel an  Avenua 
Huntington  loach,  CA  32647 

to.  H.  H.  Abramson 

Sou thus  at  lessor eh  loot l tuts 

8300  Culabra  toad 

Sam  Antonia,  TX  78284 

Dr.  ft.  C.  DeHart 

So ut hue at  ftasearch  taatitmta 

8300  Culabra  toad 

Saa  Antonio,  TX  78284 

Dr.  N.  L.  Boron 
told linger  Assoc  la tea 
110  Boat  59th  Straat 
tow  York,  Iff  10022 

to.  T.  L.  Goars 

Lockheed  Mlaal las  mad  Space  Company 
3231  tonovar  Straat 
Paid  Alto,  CA  34)04 

Mr.  Will  lam  Caywoed 
Applied  Pips tea  Laboratory 
Johna  top kins  toad 
Laurel ,  m  20110 
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Ualvaraitlaa  (Con't.) 

Profasaor  H.  W.  Liu 
Syracuse  Unlveralty 
Dapertmaat  of  Chemical  log lama  ring 
and  Metallurgy 
Syracwa,  UT  13210 

Profaaaor  8.  Bodnar 
Tochnlon  UD  Foundation 
Haifa,  Israel 

Profasaor  Warner  Goldsmith 
University  of  California 
Dapertmaat  of  Mechanical  tog laser lag 
Berks ley,  CA  94720 

Profeaaor  ft.  8.  ftlvlln 
Lehigh  Unlveralty 
Center  for  Appllcatloo 
of  Mathematics 
Bethlehem,  PA  18015 

Profeaaor  P.  A.  Cosserelll 
State  Unlveralty  of  Haw  York  at 
buffalo 

Division  of  Interdisciplinary  Studies 
Karr  Parker  Engineering  Building 
Chemistry  toed 
Buffalo,  NY  14214 

Profasaor  Joseph  L.  Rose 
Drexcl  Unlveralty 

Department  of  Mechanical  Engineering 
and  Mechanics 
Philadelphia,  PA  19104 

Professor  8.  K.  Donaldson 
University  of  Maryland 
Aaroapaca  Engineering  Department 
College  Park,  MD  20742 

Professor  Joseph  A.  Clark 
Catholic  Unlveralty  of  America 
Department  of  Mechanical  Engineering 
Washington,  D.C.  20064 
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.)  Industry  and  ftasearch  Institutes  (Coa't.) 

Or.  Robert  E.  Dunham 
Pacifica  Technology 
P.0.  Box  148 
Oel  Mar,  CA  9201 4 

Dr •  M.  P.  Kanninen 

■at telle  Coluabus  Laboratorlaa 

505  King  Avenue 

Columbus,  OH  43201 

Dr.  A.  A.  Hoc hr a  in 
Deeds la an  Associates,  loc ■ 

Spring lake  ftasearch  toad 
13110  Pradarlck  toad 
Woodbine,,  HD  21797 

Dr.  Jamas  V.  Jonas 
Swanson  Service  Corporation 
P.O.  Box  5415 

Huntington  leech,  CA  926A6 

Dr.  tobert  E.  Nlekall 
Applied  Science  sad  Technology 
3)44  North  Torroy  Pino*  Court 
Sul to  220 

La  Jolla,  CA  92037 

to.  Kovtn  Thomas 
Meetinghouse  Electric  Oorp. 

Advanced  toactora  Division 
P.O.  box  158 
Madison,  PA  15663 

Dr.  tarmard  Shaffer 
Polytechnic  Institute  of  Mew  York 
Dapt.  of  Mechanic el  aad  Aaroapaca 
toglaaerliig 
3)3  Jay  Straat 
Drooklyn,,  Iff  11021 


j 


